%% Main file 
% Calibrates the model

%% Clean up 
clear variables, close all

%% General stuff
name = 'replication';

%% Auxiliary stuff 
aux.n       = 10^3; % grid points for integration
aux.tol_G   = 10^(-10);% Tollerance for G
aux.n_m     = 10^3;% Grid for marginal product needed for inaction calculation 
aux.CN      = 0.5;% Weight in HJB iteration 

%% Exogenous parameters and moments
param.shock     = 0.01;% Size of shock (baseline has negative shock so that a positive value implies negative shock)
param.r         = log(1.05)/12;% discount rate
param.alpha     = 0.64;% degree of DRS
param.lambda    = 0.25;% job offer arrival rate
param.zeta      = 0;% exogenous separation rate
param.u_mom     = 0.06;% unemployment rate moment
param.job_to_job_mom    = 0.032;% job-to-job transition rate
param.c_mom             = 1;% hiring cost frac average monthly wage
param.lnw_gain_mom      = 0.08;% wage gain upon job-to-job transition
param.EU_cycl_mom       = -3.6;% eu cyclicality
param_ex                = param;% assign parameters directly determined

disc_list = nan(6,1);
%% loop over calibration
for i = 1:6
    
    param=param_ex;% assign exo parameters

    if i == 1
        param.inac_mom = 0;% no inaction with OJS
        parampath='parameters_OJS\';% Path for parameters
    elseif i == 2
        param.inac_mom = 0.108;% 12m employment weighted inaction rate
        parampath='parameters_10\';% Path for parameters
    elseif i == 3
        param.inac_mom = 0.168;% 12m employment weighted inaction rate
        parampath='parameters_17\';% Path for parameters
    elseif i == 4
        param.inac_mom = 0.225;% 12m employment weighted inaction rate
        parampath='parameters_22\';% Path for parameters
    elseif i == 5
        param.inac_mom = 0.271;% 12m employment weighted inaction rate
        parampath='parameters_27\';% Path for parameters
    elseif i == 6
        param.inac_mom = 0.384;% 12m employment weighted inaction rate
        parampath='parameters_38\';% Path for parameters
    end
    
    %% Initial guess
    if param.inac_mom == 0
        theta0 =[   0.15
                    0.26
                    1.0
                    0.5];
    else
        theta0 =[   0.14
                    0.5
                    1
                    0.08
                   30];
    end
    
    %% Calibrate steady state model
    theta0 = fminsearch(@ (x) moments(x,param,aux), theta0, optimset('display','iter','MaxIter',30));
    theta0 = fminsearch(@ (x) moments(x,param,aux), theta0, optimset('display','iter'));
    disc_list(i) = moments(theta0,param,aux);

    param = param_trans(theta0,param); % Assign parameters
    param = model( param, aux);% solve model given parameters
    param_ini=param;
    
    % Set bounds for cyclicality of breakdown costs
    omega_0_l=param_ini.omega_0*(1-2*param.shock);
    omega_0_u=param_ini.omega_0*1;
    
    % Calibrate cyclicality of breakdown costs
    x0 = fzero(@ (x) root_w(x,1,param,param_ini,aux), [omega_0_l omega_0_u], optimset('display','iter'));
    param_end = root_w( x0,0,param,param_ini,aux);

    %% Save parameters

    save([parampath,'param_ini_',name,'.mat'],'param_ini');
    save([parampath,'param_end_',name,'.mat'],'param_end');

    clear param param_end param_ini theta0 omega_0*

end
